set linesize 140
set rmsg on
set scheme s1mono
/******************************************************************************
	Project		:	Army Service in the All-Volunteer Era
	Author(s)	:	Kyle Greenberg	(kyle.greenberg@westpoint.edu)
					Matthew Gudgeon (matthew.gudgeon@westpoint.edu)
					Adam Isen 		(Adam.Isen@treasury.gov)
					Corbin Miller 	(Corbin.Miller@treasury.gov)
					Rich Patterson 	(rich_patterson@byu.edu)
	File Name	:	distribution.do
	Description	:	This file creates figures for the distributional outcomes.
*******************************************************************************/

*---- if not running master, set up file structure ----*
*change directory to where programs and subfolders are stored
*cd ""

if "${raw}"==""		global raw		"raw/"
if "${data}"==""	global data		"data/"
if "${output}"==""	global output	"output/"

cap mkdir	"${output}"
cap mkdir	"${output}dist/"

*---- select subprograms ----*
local regs	1
local figs	1

*---------------------*
*---- regressions ----*
*---------------------*
if `regs'==1{
	use "${data}army-treasury-analysis", clear
	
	gen birth_year=floor(firstyear-age_days)
	
	*---- create distributional outcomes ----*
	forval y = -1/19{
		local yt = `y'
		if `y' == -2 local yt m2
		if `y' == -1 local yt m1	
		di "y=`y'"
		gen tax_yr = firstyear+`y'
		
		*percentiles from a national "control" sample
		merge m:1 birth_year tax_yr male using "${raw}control_earnings", nogen keep(1 3)
		
		gen wages_0_`yt' = wages_`yt'==0 if !missing(wages_`yt')
		gen wages_0_25_`yt' = wages_`yt'>0 & wages_`yt'<=control_wage25 if !missing(wages_`yt')
		gen wages_25_50_`yt' = wages_`yt'>control_wage25 & wages_`yt'<=control_wage50 if !missing(wages_`yt')
		gen wages_50_75_`yt' = wages_`yt'>control_wage50 & wages_`yt'<=control_wage75 if !missing(wages_`yt')
		gen wages_75_100_`yt' = wages_`yt'>control_wage75 & wages_`yt'<. if !missing(wages_`yt')
		
		drop control_wage* tax_yr
	}
	
	keep firstafqt firstafqt_fy firstyear inst31 k31* instk31* inst50 k50* instk50* access quarterFE wages_0_*  wages_0_25_*  wages_25_50_* wages_50_75_* wages_75_100_*
	
	*snapshot 1
	snapshot save
	
	foreach outcome in wages_0 wages_0_25 wages_25_50 wages_50_75 wages_75_100{
		*---- initialize output file ----*
		clear
		set obs 22
		gen y = _n - 3
		foreach cut in 31 50{
			foreach x in b se N m{
				gen `x'_`outcome'_`cut' = .
			}
		}
		save "${output}dist/dist_results_`outcome'", replace
		
		forval y = -1/19{
			local yt = `y'
			if `y' == -2 local yt m2
			if `y' == -1 local yt m1	
			di "y=`y'"
			*---- run regressions ----*
			snapshot restore 1
			
			*1st cutoff (AFQT=31)
			ivregress 2sls `outcome'_`yt' k31 instk31 k31_2 instk31_2 (access=inst31) ib175.quarter if inrange(firstafqt,12,49), robust
			local b_31 = _b[access]
			local se_31 = _se[access]
			local N_31 = e(N)
			sum `outcome'_`yt' if e(sample)
			local m_31 = r(mean)
				
			*2nd cutoff (AFQT=50)
			ivregress 2sls `outcome'_`yt' k50 instk50 k50_2 instk50_2 (access=inst50) ib175.quarter if inrange(firstafqt,31,68), robust
			local b_50 = _b[access]
			local se_50 = _se[access]
			local N_50 = e(N)
			sum `outcome'_`yt' if e(sample)		
			local m_50 = r(mean)
			
			*---- write estimates to dataset ----*
			use "${output}dist/dist_results_`outcome'", clear
			foreach cut in 31 50{
				foreach x in b se N m{
					replace `x'_`outcome'_`cut' = ``x'_`cut'' if y == `y'
				}
			}
			save "${output}dist/dist_results_`outcome'", replace
		}
	}
	snapshot erase 1
}	//end if
else di "*---- skipping distributional regressions ----*"

*------------------------*
*---- output figures ----*
*------------------------*
if `figs'==1{
	use "${output}dist/dist_results_wages_0", clear
	gen hi_wages_0_31 = b_wages_0_31 + 1.96*se_wages_0_31
	gen lo_wages_0_31 = b_wages_0_31 - 1.96*se_wages_0_31
	gen hi_wages_0_50 = b_wages_0_50 + 1.96*se_wages_0_50
	gen lo_wages_0_50 = b_wages_0_50 - 1.96*se_wages_0_50
	foreach outcome in wages_0_25 wages_25_50 wages_50_75 wages_75_100 {
		merge 1:1 y using "${output}dist/dist_results_`outcome'", nogen
		gen hi_`outcome'_31 = b_`outcome'_31 + 1.96*se_`outcome'_31
		gen lo_`outcome'_31 = b_`outcome'_31 - 1.96*se_`outcome'_31
		gen hi_`outcome'_50 = b_`outcome'_50 + 1.96*se_`outcome'_50
		gen lo_`outcome'_50 = b_`outcome'_50 - 1.96*se_`outcome'_50	
	}
	
	foreach cut in 31 50{
		*just wages_0
		tw (line b_wages_0_`cut' y, lcolor(black) lwidth(medthick))	///
			(rcap hi_wages_0_`cut' lo_wages_0_`cut' y, lcolor(black) lp(longdash) lw(thin)),	///
			xline(0, lpattern(dash) lcolor(gs8)) yline(0, lpattern(dash) lcolor(gs8))	///
			ytitle("Probability of Earnings w/in Given Range") ytick(-.2(.1).4) ylabel(-0.2(0.1)0.4, angle(0) format(%4.2f))	///
			xtitle("Years Since Application") xtick(-1(1)19) xlabel(-1(2)19)	///
			legend(order(1 "Zero Wages"))
		graph export "${output}dist/dist_wages_0_`cut'.pdf", as(pdf) replace

		*add wages_0_25
		tw (line b_wages_0_25_`cut' y, lcolor(black) lwidth(medthick))	///
			(rcap hi_wages_0_25_`cut' lo_wages_0_25_`cut' y, lcolor(black) lp(longdash) lw(thin))	///
			(line b_wages_0_`cut' y, lcolor(gs4) lwidth(medthick) lpattern(dot)) , ///
			xline(0, lpattern(dash) lcolor(gs8)) yline(0, lpattern(dash) lcolor(gs8))	///
			ytitle("Probability of Earnings w/in Given Range") ytick(-.2(.1).4) ylabel(-.2(0.1).4, angle(0) format(%4.2f))	///
			xtitle("Years Since Application") xtick(-1(1)19) xlabel(-1(2)19)	///
			legend(order(3 "Zero Wages" 1 "Wages 0-25 percentile"))
		graph export "${output}dist/dist_wages_0_25_`cut'.pdf", as(pdf) replace

		*add wages_25_50
		tw (line b_wages_25_50_`cut' y, lcolor(black) lwidth(medthick))	///
			(rcap hi_wages_25_50_`cut' lo_wages_25_50_`cut' y, lcolor(black) lp(longdash) lw(thin))	///
			(line b_wages_0_25_`cut' y, lcolor(maroon) lwidth(medthick) lpattern(dot)) ///
			(line b_wages_0_`cut' y, lcolor(gs6) lwidth(medthick) lpattern(dot)), ///
			xline(0, lpattern(dash) lcolor(gs8)) yline(0, lpattern(dash) lcolor(gs8))	///
			ytitle("Probability of Earnings w/in Given Range") ytick(-.2(.1).4) ylabel(-.2(0.1).4, angle(0) format(%4.2f))	///
			xtitle("Years Since Application") xtick(-1(1)19) xlabel(-1(2)19)	///
			legend(order( 4 "Zero Wages" 3 "Wages 0-25 percentile"  1 "Wages 25-50 percentile" ))
		graph export "${output}dist/dist_wages_25_50_`cut'.pdf", as(pdf) replace

		*add wages_50_75
		tw (line b_wages_50_75_`cut' y, lcolor(black) lwidth(medthick))	///
			(rcap hi_wages_50_75_`cut' lo_wages_50_75_`cut' y, lcolor(black) lp(longdash) lw(thin))	///
			(line b_wages_25_50_`cut' y, lcolor(red) lwidth(medthick) lpattern(dot)) ///
			(line b_wages_0_25_`cut' y, lcolor(maroon) lwidth(medthick) lpattern(dot))	///
			(line b_wages_0_`cut' y, lcolor(gs6) lwidth(medthick) lpattern(dot)), ///
			xline(0, lpattern(dash) lcolor(gs8)) yline(0, lpattern(dash) lcolor(gs8))	///
			ytitle("Probability of Earnings w/in Given Range") ytick(-.2(.1).4) ylabel(-.2(0.1).4, angle(0) format(%4.2f))	///
			xtitle("Years Since Application") xtick(-1(1)19) xlabel(-1(2)19)	///
			legend(order(5 "Zero Wages" 4 "Wages 0-25 percentile"  3 "Wages 25-50 percentile"  1 "Wages 50-75 percentile"))
		graph export "${output}dist/dist_wages_50_75_`cut'.pdf", as(pdf) replace

		*add wages_75_100
		tw (line b_wages_75_100_`cut' y, lcolor(black) lwidth(medthick))	///
			(rcap hi_wages_75_100_`cut' lo_wages_75_100_`cut' y, lcolor(black) lp(longdash) lw(thin))	///
			(line b_wages_50_75_`cut' y, lcolor(orange) lwidth(medthick) lpattern(dot)) ///
			(line b_wages_25_50_`cut' y, lcolor(red) lwidth(medthick) lpattern(dot)) ///
			(line b_wages_0_25_`cut' y, lcolor(maroon) lwidth(medthick) lpattern(dot))	///
			(line b_wages_0_`cut' y, lcolor(gs6) lwidth(medthick) lpattern(dot)), ///
			xline(0, lpattern(dash) lcolor(gs8)) yline(0, lpattern(dash) lcolor(gs8))	///
			ytitle("Probability of Earnings w/in Given Range") ytick(-.2(.1).4) ylabel(-.2(0.1).4, angle(0) format(%4.2f))	///
			xtitle("Years Since Application") xtick(-1(1)19) xlabel(-1(2)19)	///
			legend(order(6 "Zero Wages" 5 "Wages 0-25 percentile"    4 "Wages 25-50 percentile"  3 "Wages 50-75 percentile" 1 "Wages 75-100 percentile"))
		graph export "${output}dist/dist_wages_75_100_`cut'.pdf", as(pdf) replace

		*all at once, no SE bars
		tw (line b_wages_75_100_`cut' y, lcolor(black) lwidth(medthick) lp(dash))	///
			(line b_wages_50_75_`cut' y, lcolor(blue) lwidth(medthick) lpattern(dash)) ///
			(line b_wages_25_50_`cut' y, lcolor(orange) lwidth(medthick) lpattern(dot))	///
			(line b_wages_0_25_`cut' y, lcolor(maroon) lwidth(medthick) lpattern(dot)) ///
			(line b_wages_0_`cut' y, lcolor(gs6) lwidth(medthick) lpattern(dot)), ///
			xline(0, lpattern(dash) lcolor(gs8)) yline(0, lpattern(dash) lcolor(gs8))  xtick(-1(1)19) xlabel(-1(2)19) legend( order( 5 "Zero Wages" 4 "Wages 0-25 percentile" 3 "Wages 25-50 percentile" 2 "Wages 50-75 percentile" 1 "Wages 75-100 percentile"))
		graph export "${output}dist/dist_all_`cut'.pdf", as(pdf) replace
	}
}	//end if
else di "*---- skip outputting distributional figures ----*"

